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Due to a large discrepancy between theory and experiment, the electronic character of crystalline 
boron carbide B13C2 has been a controversial topic in the field of icosahedral boron-rich solids. We 
demonstrate that this discrepancy is removed when configurational disorder is accurately considered 
in the theoretical calculations. We find that while ordered ground state B13C2 is metallic, config¬ 
urationally disordered B13C2, modeled with a superatom-special quasirandom structure method, 
goes through a metal to non-metal transition as the degree of disorder is increased with increasing 
temperature. Specifically, one of the chain-end carbon atoms in the CBC chains substitutes a neigh¬ 
boring equatorial boron atom in a B12 icosahedron bonded to it, giving rise to a Bn(BBC) unit. 

The atomic configuration of the substitutionally disordered B13C2 thus tends to be dominated by 
a mixture between Bi2(CBC) and BiiC^(BBC). Due to splitting of valence states in BiiC^(BBC), 
the electron deficiency in Bi2(CBC) is gradually compensated. 


I. INTRODUCTION 

Boron carbide, a covalent solid, is a promising candi¬ 
date for a wide range of applications, e.g. lightweight 
armor plates, cutting tools, thermoelectric conversion 
and a neutron absorbing material in nuclear reactors [1]. 
Boron carbide has also been realized in a neutron detec¬ 
tor application due to a high cross section for thermal 
neutron reaction of ^^B [2-5]. Experimental findings [6- 
9] reveal a fundamental crystalline structure of boron 
carbide with 12-atom icosahedra placed at vertices of 
a rhombohedral unitcell {RZm space group) coexisting 
with a 3-atom intericosahedral chain aligning itself along 
the longest diagonal in the rhombohedron. Except 
from a small shift in lattice spacing, the crystalline 
structure of boron carbide is unchanged over a wide 
range of compositions (~8-20 at.% C), corresponding to 
a single-phase region [1, 10]. Thus, a configurationally 
disordered solid solution of boron and carbon atoms 
within the structural units are inevitable to take place. 
Consequently, to get an understanding of how such a 
substitutional disorder influences the properties of boron 
carbide, knowing how boron and carbon atoms are 
distributed to form boron carbide is of importance for its 
technological applications. Unfortunately, experimen¬ 
tally identifying the exact atomic positions of carbon 
atoms in these compounds, e.g. by x-ray diffraction 
technique, is extremely difficult due to the fact that 
the atomic form factors of boron and carbon are very 
similar. This has become a controversial issue about 
the structural changes of boron carbide as the carbon 
concentration changes, in particular for B13C2 with 
-13.33 at.% C [11]. 

It is generally accepted that at the carbon-rich 
limit, i.e. B4C with 20 at.% C, the structural units of 


boron carbide are dominated by Bn(CBC), where 
p denotes the polar sites of icosahedra [12-14]. This 
presumption is also in line with results from our recent 
study of configurational disorder in B4C [15]. As the 
carbon concentration decreases, approaching B13C2, 
the replacement of carbon by boron can take place 
either within the icosahedra or in the chains, resulting 
in two different models: Bi2(CBC) or BiiC(BBC), 
respectively, for B13C2. The former model is evident 
mainly by ab initio calculations [16-18], where it is 
shown to have considerably lower energy as compared 
to the BiiC(BBC) [16, 19]. Meanwhile the latter model 
is consistent with the analyses of structural data from 
x-ray diffraction [20] and Raman spectra [21-23] of 
boron carbide at different at.% C in its single-phase 
region. 

Moreover, the calculations predict that the presumed 
Bi 2(CBC) is metallic [16, 24, 25] due to its electron 
deficiency (1 hole/unitcell). The experimental obser¬ 
vations [26, 27] however reveal that, throughout the 
single-phase region, boron carbide is a semiconductor, 
thus giving rise to a discrepancy between theory and 
experiment. According to a study of the bonding 
nature of B13C2 by Balakrishnarajan et al [28], they 
suggested that the semiconducting behavior in B13C2 
originates from unavoidable boron/carbon substitutional 
disorder, leading to localization of electronic states. 
This seems corresponding to the statement proposed 
by Schmechel [29] and Werheit [30] that the electron 
deficiency in B13C2 can be compensated by defects, 
e.g., substitutional defects and vacancies, splitting off 
some valence states into a band-gap. However, no 
detailed suggestion or calculations demonstrating how 
such effects should be realized were present. A recent 
theoretical study conducted by Shirai et al [19] proposed 
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structural models for B 13 C 2 , consisting of BiiC^(CBC), 
Bi 2 (CBC), and Bi 2 (B 4 ) units. According to their 
models, the presence of the Bi 2 (B 4 ) units can partially 
solve the band-gap problem and the hole density reduces 
to 0.25 hole/unitcell, indicating less metallic character. 
However, their models did not demonstrate the semi¬ 
conducting character for B 13 C 2 , the effects of specific 
configurations between their suggested structural units 
were not fully considered, and no thermodynamic sta¬ 
bility analysis including entropic effects were performed. 
Accordingly, further investigations not only to resolve 
the inconsistencies between experimental and theoretical 
studies of B 13 C 2 , especially the electronic structure, but 
also to understand the influence of substitutional disor¬ 
der on the properties of boron carbide, are necessary. 

In this work, we investigate the phase stability, 
within a mean-field approximation for the configura¬ 
tional entropy, and the properties of substitutionally 
disordered B 13 C 2 using first-principles calculations. A 
random alloy theory-based superatom-special quasir¬ 
andom structure (SA-SQS) method [15, 31] is used to 
model disordered configurations of high concentrations 
of low-energy defects in B 13 C 2 , as previously applied 
to B 4 C [15]. We demonstrate that, at elevated tem¬ 
perature, B 13 C 2 is thermodynamically favorable to 
substitutionally disorder so that the Bi 2 (CBC) units 
coexist with the BnC®(BBC) units, where e stands for 
equatorial sites. We then show that the existence of the 
BnC®(BBC) units compensates the electron deficiency 
in the Bi 2 (CBC) units. As the degree of substitutional 
disorder increases with temperature under equilibrium 
condition, a transition from the ordered Bi 2 (CBC) 
metallic state to the disordered B 13 C 2 semiconducting 
state is achieved. 


II. COMPUTATIONAL DETAILS 

All calculations are performed within a realm of Den¬ 
sity functional theory (DFT) as implemented in the 
Vienna ab initio simulation package (VASP) [32, 33]. 
For the total-energy calculations, the projector aug¬ 
mented wave (PAW) method [34] with the Perdew-Becke- 
Ernzerhof (PBE96) generalized gradient approximation 
(GGA) [35], as the exchange-correlation functional, is 
used. Regarding the equilibrium volume optimization, 
the internal degrees of freedom and the cell shape are 
fully relaxed for a set of different fixed volume calcu¬ 
lations. Eor the density of states calculations, we em¬ 
ploy three different exchange-correlation functionals, i.e. 
GGA-PBE96, MBJ-GGA [36, 37], and HSE06 [38, 39]. 
Energy cutoff of 400 eV is used for the plane-wave cal¬ 
culations. Meanwhile for the Brillouin zone integration, 
Monkhorst-Pack k-point mesh [40] is chosen. Since var¬ 
ious supercell sizes are used in this work, in order to 


TABLE I: Different k-point grids for the Brillouin zone inte¬ 
gration, which are used m:(l) Cell optimization calculations 
and (2) Density of states calculations, for different supercell 
sizes. 


Cell size (Aatoms) 

k-point grids 

Cell optimizations 

Density of states 

1x1x1 (15) 

5x5x5 

9x9x9 



5x5x5 (HSE06) 

2x1x1 (30) 

5x5x5 

9x9x9 

2x2x2 (120) 

5x5x5 

9x9x9 

4x3x3 (540) 

3x3x3 

5x5x5 

4x4x3 (720) 

3x3x3 

3x3x3 


obtain a good accuracy within a reasonable computa¬ 
tional time for different kinds of calculations, different 
k point grids are used and given in TABLE 1. Visual¬ 
izations of atomic configurations are obtained with the 
VESTA package [41]. 


III. DILUTE SUBSTITUTIONAL DEFECTS 

Due to the complexity of the 15-atom structural unit 
and the similarity of boron and carbon atoms, numer¬ 
ous kinds of substitutional defects are conceivable. Tak¬ 
ing them all into account in modelling a substitutionally 
disordered B 13 C 2 is a difficult task. It is thus neces¬ 
sary to single out defects, which are low in energy and 
likely to appear in the structure at high concentrations. 
In this way, substitutionally disordered configurations 
of B 13 C 2 can be properly modelled. We first consider 
different kinds of dilute substitutional defects, involved 
with only 2 or 4 atoms, in a 2x2x2 supercell consist¬ 
ing of eight Bi 2 (CBC) units, which is believed to be the 
configurational ground state of B 13 C 2 (see Fig. 1 ). We 
note that substitutional defects, considered in this sec¬ 
tion, are generated only by swapping boron and carbon 
atoms within structural units, which does not alter the 
global stoichiometry. Thus, all defective structures have 
the stoichiometry of B 13 C 2 . Substitutional defect forma¬ 
tion energies, denoted by AE^^e/ect, are calculated by the 
equation; 

^^defect — ^defect ^Bi2{C BC) •> (f) 

where ^defect and ^Bi2{cbc) are defined as the total en¬ 
ergy of the defective structure and the ground state to¬ 
tal energy, respectively. Eor each defective structure, we 
find that volume optimization of the supercell reduces 
/SFjdefect by less than 5%. The volume is thus kept fixed 
at the equilibrium volume of the ground state, and during 
the total energy calculation, only the internal degrees of 
freedom are allowed to relax. AE^^g/ect for different types 
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of substitutional defects are listed in TABLE 11. We note 
that all considered defective structures are less energeti¬ 
cally favorable as compared to Bi 2 (CBC), but still, in the 
considered dilute concentration, stable against a phase 
decomposition into a-boron and diamond. 



FIG. 1: A 12 -atom icosahedron and two 3-atom intericosahe- 
dral chains, illustrating the structural units of the predicted 
ground state of B 13 C 2 . Light and dark spheres represent 
boron and carbon atoms, respectively. The notations p{l- 
6 ) and e(l- 6 ) stand for two different crystallographic sites in 
the icosahedron, i.e. polar (bonded to neighboring icosahe- 
dra) and equatorial (bonded to intericosahedral chains) sites, 
respectively. On the one hand, the notations c(l-3) denote 
the chain sites, which are, in fact, equivalent to c(4-6). 

The atomic positions of carbon and boron atoms, in 
this work, are denoted by superscripts p, e, and c, cor¬ 
responding to their positions at polar, equatorial, and 
chain sites, respectively, and the numbers denote the po¬ 
sitions according to the notations in Fig. 1. For example, 
B(c 2 ) refers to a boron atom residing at the center po¬ 
sition of the chain unit, i.e. the c2 position in Fig. 1. 
From TABLE II, we see that a swap of by B^®^^ 

yields a substitutional defect with the lowest AEdefect 
of 0.605 eV, i.e. a substitution of a chain-end carbon 
atom (cl) for a neighboring equatorial boron atom (el) 
in the icosahedron, denoted by BnC®(BBC). This for¬ 
mation energy is increased by approximately a factor of 
1.5, if the substitution of the other chain-end carbon 
atom (c3) for a chain-center boron atom (c2) is taken 
into account, i.e. BiiC®(BCB) with AEdefect = 0.911 
eV. Instead of a substitution of for as in the 

case of BnC®(BBC), AEdefect slightly further increases 
as compared to that of BnC®(BCB), if the substitution 
of a chain-end carbon atom occurs at a polar site in 
the neighboring icosahedron, for example at the pi po¬ 
sition {AEdefect = 0.935 eV for a substitution of C^^^^ 


TABLE II: Substitutional defect formation energies, 
AEdefect: ^ 2 x 2 x 2 supcrcell with respect to the non¬ 

defective structure of B 13 C 2 , i.e. Bi 2 (CBC). The super¬ 
scripts of the notations in the parentheses correspond to the 
positions given in FIG. 1. 

Defective structure 

AEdefect (eV/defect) 

Non-defective structure 


Bi 2 +(C-B-C) 

0 

Disordered chain 


Bi 2 +(B('=i^-C('= 2 )-C) 

2.570 

Bi2+(C-C('=2^-B('=3)) 

2.570 

Polar sites 


BiiC(p^)+(B(‘=^)-B-C) 

1.173 

BiiC(p^)+(C-B-B(=^)) 

0.935 


3.812 

Equatorial sites 


BiiC(®^>+(B(‘=i>-B-C) 

0.605 

BiiC(®^)+(C-B-B('='^)) 

1.548 

BiiC(®®^+(B('=i)-B-C) 

1.666 


0.911 

BiiC(®^)+BiiC(®^)+(B(‘=^>-B-B(‘='^)) 

3.530 

Bipolar defect 



1.517 

BioC 2 ^^^’^^^ + (B('= 1 ’'=®)-B-C )2 

2.027 

BioC 2 ^P’^’P®^(B(‘= 1 ’=«>-C-C )2 

2.185 

Biequatorial defect 


BioC2<'^^’"^4(B(=1’‘=®>-C-C)2 

1.585 


for B^^^\ yielding BnC^(CBB) and AEdefect = 1-173 eV 
a substitution of C^^^^ for yielding BnC^(BBC)). 

On the other hand, a substitution of for B^®^\ or 

BnC®(CBB), yields even higher AEdefect^ he. 1.548 eV, 
because of, in this case, a formation of a C-C bond be¬ 
tween the chain unit and the icosahedron. This indicates 
that a direct chemical bonding between carbon atoms is 
unfavorable in this compound in line with the earlier find¬ 
ing for B4C [15]. This is also supported by a substitution 
of C^^^^ for B^®^^ {AEdefect = 1.666 eV), resulting in the 
same kind of C-C bond as in the case of BiiC®(CBB). 
Another support is a very high AEdefect of 2.570 eV in 
the case of disordered chain, i.e. Bi 2 (BCC), where there 
exists an intrachain C-C bond. Since the BnC®(BBC) 
type of substitutional defect is found to have the lowest 
AEdefect: with Other types of substitutional defects hav¬ 
ing considerably higher AEdefect: f^is type of defect is 
likely to dominate disordered configurations of B 13 C 2 . 


IV. SPLITTING OF VALENCE STATES 

We calculate the electronic density of state (DOS) for 
a 15-atom unitcell of Bi 2 (CBC), as shown in Eig. 2(a). 
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FIG. 2: (Color online) Density of states of (a) Bi 2 (CBC), 
(b) BiiC^(BBC), (c) Bi 2 (BBC), and (d) a combination of 
Bi 2 (CBC) and Bn(BBC). The black solid line, the red 
dashed line, and the blue dashed-dotted line indicate the 
electronic DOS, obtained by using the GGA-PBE96, the 
MBJ-GGA, and the HSE06, respectively, for the exchange- 
correlation functional. The highest occupied state is located 
at 0 eV for all cases. 


The result, obtained from the GGA-PBE96 functional, 
shows the Fermi level is located in the valence band, thus 
resulting in a metallic behavior of the Bi 2 (CBC) unit, 
due to its electron deficiency. This corresponds to the 
theoretical results, reported in the literature [16, 24, 25]. 
However, experimentally, boron carbide is a semiconduc¬ 
tor throughout the single-phase region [26, 27]. We then 
use the MBJ-GGA and the HSE06 functionals, which 
are known to give a better description for electronic 
DOS, in particular the band-gap of semiconductors, 
compared to the conventional GGA-PBE96. The DOS, 
obtained from the latter two functionals, confirm that 
the Bi 2 (CBC) unit are metallic. We thus conclude that 
the metallic character of the Bi 2 (CBC) unit is, in this 
study, independent of the exchange-correlation function¬ 
als and the discrepancy in electronic properties between 
experiment and theory does not arise from inaccuracies 
in the used exchange-correlation functionals. 

Besides, we calculate the DOS of BnC®(BBC) in a 
15-atom unitcell, as shown in Fig. 2(b). In this case, we 
observe splitting of two valence states into a band-gap 
and the Fermi level is located in the midpoint of the 
splitting states. Due to this valence band splitting, we 
consider also spin polarization in the Bi 2 (CBC) and the 
BnC®(BBC) units, using the HSE06 functional. The 
results show non-zero magnetic moment in both cases. 
The total energy of the spin-polarized BnC®(BBC) is 
lower by 0.3 eV/unitcell, as compared to that calculated 
in the case of non-spin polarized calculation. On the 
contrary, the non-magnetic solution gives the lowest 
total energy to the Bi 2 (CBC) unit. These results 
indicate that it might be necessary to take into account 


the effect of spin polarization, when calculating boron 
carbide configurations where the Fermi level falls into 
narrow defect-like states. 

Unlike the BnC®(BBC) unit, the DOS of the 
BnC®(BCB) unit (not shown) shows no splitting of 
valence states and the Fermi level lies within the 
valence band as in the case of the Bi2(CBC) unit. 
Similarly to the Bi 2 (CBC) unit, the BiiC®(BBC) and 
the BiiC®(BCB) units are both metallic. We suggest 
that such a valence band splitting in the BnC®(BBC) 
unit originates from the presence of the BBC chain. Our 
suggestion is supported by the valence band splitting in 
the boron-rich Bi 2 (BBC) unit, as shown by Fig. 2(c). 
It is noted that the Bi 2 (BBC) unit is a semiconductor, 
since the Bi 2 (BBC) unit has one less electron, compared 
to the BnC®(BBC) unit, thus resulting in the two empty 
splitting valence states in the Bi 2 (BBC) unit. The same 
kind of valence band splitting is also observed in the 
BnC^(CBB), the BnC^(BBC), and the BnC^(CBB) 
units (not shown). 

We further investigate the electronic properties 
of a combination between two different units, i.e. 
Bi 2 (CBC) and BnC®(BBC), in a 2 x 1 x 1 supercell (30 
atoms). As illustrated by Fig. 2(d), the supercell of 
B 13 C 2 exhibits a semiconducting character with a 
band-gap of 1.36 (2.09) eV according to a use of the 
GGA-PBE96 (the MBJ-GGA) functional. Also the 
spin polarized calculation, using the HSE06 functional 
reveals that it is non-magnetic. This is because an 
electron in the splitting states of the BnC®(BBC) unit 
occupies an empty state in the Bi 2 (CBC) unit, thus 
compensating the electron deficiency. This observation 
is in line with the speculations by Balakrishnarajan [28], 
Schmechel [29], and Werheit [30] that boron/carbon 
substitutional disorder could lead to localization of 
electronic states, compensating electron deficiency in 
B 13 C 2 . Such a compensation of electron deficiency 
can also take place in a combination of Bi 2 (CBC) and 
BiiC^(CBB) (BiiC^(BBC)) with a band-gap of 1.0 
( 1 . 8 ) eV (not shown). Since neither the Bi 2 (CBC) 
unit nor its combination with the BnC®(BBC) unit is 
magnetic, the results presented in the rest of this work 
will be based only on the non-spin polarized calculations. 


V. MODELING CONFIGURATIONALLY 
DISORDERED B13C2 

We demonstrated in section III that the formation en¬ 
ergies of the BiiC^(BCB), BiiC^(CBB), BiiC^(BBC), 
and BiiC®(CBB) units in a matrix of Bi2(CBC) are 
somewhat larger than that of the BnC®(BBC) unit 
(see TABLE II). In section IV, we also showed that the 
electron deficiency in B 13 C 2 can be compensated by a 
mixture of Bi2(CBC) and BnC^(BBC) with a ratio of 
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FIG. 3: A superatom basis for modelling configurationally disordered B 13 C 2 . Notations and colors indicate, respectively, the 
same crystallographic sites and the same kinds of atoms as described in Fig. 1. The types of superatoms depend on the positions 
of the two carbon atoms residing in that superatom. A superatom (a) stands for the Bi 2 (CBC) unit, while the Bn(BBC) 
unit can be represented, for example, by a superatom (b) or (c). 


1 : 1 , keeping the B 13 C 2 stoichiometry. Consequently, 
snbstitntionally disordered configurations of B 13 C 2 are 
modelled here by taking these findings into consideration. 
The superatom-special quasirandom structure (SA-SQS) 
approach [15] is used to construct configurationally 
disordered B 13 C 2 with high concentrations of low energy 
defects distributed in a manner following the SQS 
approach for alloy theory [31]. Fig. 3 represents a chosen 
superatom basis that is used to model snbstitntionally 
disordered B 13 C 2 , mainly based on a combination of 
Bi 2 (CBC) and BnC®(BBC). In this case, the superatom 
is focused on the chain, in which the equatorial sites 
belonging to a single icosahedron, within the basis, 
are replaced by the corresponding equatorial sites (e 2 
to e 6 in Fig. 3) of neighboring icosahedra bonded to 
the original intericosahedral chain. This superatom 
basis allows for constructions of both the Bi 2 (CBC) 
ground state unit, as well as all six displacements of the 
chain-end carbon atoms to their neighboring equatorial 
sites found to be low in AEdefect^ i-e. the BnC®(BBC) 
unit. 

A superatom, illustrated in Fig. 3(a), stands for 
the Bi 2 (CBC) unit, meanwhile to define a superatom 
representing the BiiC®(BBC) unit, we assign some 
constraints to the two carbon atoms within the basis 
in order to avoid the C-C bond between the chain unit 
and the icosahedron within a superatom. That is the 
chain-end carbon atoms can swap their positions only 
with the neighboring equatorial boron atoms bonded to 
them and only one of them is allowed to swap within the 
single superatom. That is, can swap its position 

either with B^® 1 )^ B^®^\ or B^®^^ (Fig. 3(b)), meanwhile 
C^^^^ is allowed to swap its position with B^®^\ or 

(Fig. 3(c)). Consequently, there can be as many 
as seven types of superatoms considered in modelling 
configurationally disordered B 13 C 2 , consisting of the 


Bi 2 (CBC) and the BiiC®(BBC) units. 

In this study, different snbstitntionally disordered 
configurations of B 13 C 2 are considered within 4x3x3 
(540 atoms) and 4x4x3 (720 atoms) supercells. Results 
from the calculations of some selected configurations 
are shown in TABLE III. The notation for each con¬ 
figuration, e.g., e(l), e(l, 4), etc., refers to the atomic 
positions within the icosahedra, corresponding to those 
given in Fig. 1, which are substituted by carbon atoms 
from the intericosahedral chains. We note that, to com¬ 
pensate the electron deficiency, a concentration of the 
Bi 2 (CBC) unit is fixed at 50% in all cases. Meanwhile 
the other 50% is attributed to those with the BBC (or 
CBB) chain units, e.g., BiiC^(BBC), BiiC^(BBC), and 
BiiC^(CBB). The latter 50% can be equally divided, 
according to types of superatoms representing the units 
with the BBC (or CBB) chain. For example, in the e(l- 
3) configuration, ~16.67% of equatorial boron atoms, at 
el, e2, and e3 positions, are equally substituted in a 
quasirandom manner by their neighboring carbon atoms 
at the cl position. In this case, there are three types of 
the BiiC®(BBC) unit, i.e. ~16.67% for each. We note 
that in the e(l- 6 )+p(l- 6 ) configuration, we define twelve 
more types of superatoms, respecting to substitution 
of the chain-end carbon atom (C^"^^^ or C(^^)) for a 
polar boron atom within the single superatom. This 
is to allow not only the substitution of the chain-end 
carbon atoms for the equatorial boron atoms, but also 
for the polar boron atoms, i.e. to take into account 
the BiiC^(CBB) and the BiiC^(BBC) units, found 
in section III to be reasonably low in AEdefect and of 
relevance for the carbon-richer compositions approach¬ 
ing B 4 C [15], to model snbstitntionally disordered B 13 C 2 . 
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TABLE III: Energy difference, formation energy with respect to a-boron and diamond, supercell size, lattice constant, and 
electronic band-gap of different substitutionally disordered configurations of B 13 C 2 . In the first column, Bi 2 (CBC), as a 
reference, stands for the ordered-ground state of B 13 C 2 . The notations e and p stand for the substitutionally disordered B 13 C 2 
with a replacement of boron atoms in the icosahedra by carbon atoms in the intericosahedral chains, respectively, at different 
equatorial and polar sites, indicated by the numbers, corresponding to the positions in Eig. 1, in the parenthesis. AE, in the 
second column, denotes the energy difference per superatom (s.a.) with respect to the ordered state. in the third 

column, denotes the formation energy, calculated from a-boron and diamond. The lattice parameter a and the angle a are 
rhombohedrally averaged. E^ denotes the electronic band-gap obtained from GGA-PBE96 (MBJ-GGA). 


Gonfiguration 

AE 

(eV/s.a.) 

Ae/°™ 
(eV / atom) 

Supercell size 
(Number of atoms) 

a 

(A/unitcell) 

a 

n 

E. 

(eV) 

Bi2(CBC) 

0 

-0.092 

1x1x1 (15) 

5.199 

65.90 

metal (metal) 

Bi2(CBC) 

0 

-0.092 

4x3x3 (540) 

5.197 

65.98 

metal 

Bi2(CBC) “ 

0 

-0.094 

1x1x1 (15) 

5.196 

65.54 

metal 

BiiC"(BBC) 

1.778 

0.026 

1x1x1 (15) 

5.166 

66.20 

metal (metal) 

BiiC"(BBC) “ 

2.093 

0.045 

1x1x1 (15) 

5.167 

66.02 

- 

e(l) 

0.256 

-0.075 

4x3x3 (540) 

5.169 

65.85 

1.04 

e(l, 4) 

0.248 

-0.076 

4x3x3 (540) 

5.167 

65.84 

0.84 

e(l-3) 

0.253 

-0.076 

4x3x3 (540) 

5.169 

65.84 

1.16 

e(l-6) 

0.258 

-0.075 

4x3x3 (540) 

5.169 

65.85 

1.12 (1.65) 

e(l-6)+p(l-6) 

0.420 

-0.065 

4x4x3 (720) 

5.159 

66.08 

0.88 

Exp. * 

- 

- 

- 

5.196 

65.62 

- 

Exp. 

- 

- 

- 

5.185 

65.60 

- 

Exp. 

- 

- 

- 

- 

- 

2.09 


®Ref. [16] - Bylander et al. (LDA-PP) 

^Ref. [42] - Kirfel et al. (XRD) 

^Ref. [8] - Morosin et al. (XRD) 

^Ref. [43] - Werheit et al. (Photoluminescence) 


VI. RESULTS AND DISCUSSION 


a. Properties of disordered B13C2 


the experimental observations that B 13 C 2 is a semicon¬ 
ductor. 


From TABLE III, the BnC®(BBC) unit has consider¬ 
ably higher energy as compared to the Bi 2 (CBC) unit. 
This is in good agreement with the previous calcula¬ 
tion [16]. Except the BiiC®(BBC) unit, the Bi 2 (CBC) 
unit and all substitutionally disordered configurations of 
B 13 C 2 are stable against a phase decomposition into a- 
boron and diamond. Interestingly, all of the disordered 
B 13 C 2 reveal a semiconducting character with a GGA- 
PBE96-based electronic band-gap approximately rang¬ 
ing between 0.8 and 1.2 eV. By taking into account the 
fact that the GGA-PBE96 functional regularly underes¬ 
timates real band-gaps, our results are in line with the 
experimental band-gap of 2.09 eV, measured by photo¬ 
luminescence [43]. A better agreement of the calculated 
band-gap with the experiment can be achieved by a use of 
other exchange-correlation functionals that give a better 
description of electronic DOS. For instance, the band- 
gap of the e(l- 6 ) configuration becomes 1.65 eV with a 
use of MBJ-GGA. The electronic DOS of the e(l- 6 ) and 
e(l- 6 )+p(l- 6 ) configurations are shown in Fig. 4. It is 
worth noting that the presented results are different to 
all previous calculations for B 13 C 2 , and correspond to 



FIG. 4: (Golor online) Density of states of (a) the e(l-6), and 
(b) the e(l-6)+p(l-6) configurations. The black solid line, 
and the red dashed line, indicate the electronic DOS, obtained 
by using the GGA-PBE96, and the MBJ-GGA, respectively 
for the exchange-correlation functional. The highest occupied 
state is located at 0 eV. 

Experimentally, boron carbide has the rhombohedral 
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symmetry {R3m). The lattice parameters (a, 6 , c) and 
also the angles (a, /3, 7 ) at equilibrium for all disordered 
configurations listed in TABLE III are slightly deviating 
from each other by less than 1 %. We attribute this devi¬ 
ation to substitutional disorder within the supercell with 
a finite size. The existence of the icosahedral carbon in 
the BiiC®(BBC) unit results in the monoclinic distor¬ 
tion as does the ground state B 4 C, i.e. BiiC^(CBC). As 
a result, the lattice parameters and the angles, given in 
TABLE III are rhombohedrally averaged. The averaged 
parameters a and a for disordered B 13 C 2 are slightly un¬ 
derestimated by less than 0.5% and 0.3%, respectively, 
when compared to the experiments [ 8 , 42]. The averaged 
a and a for the BnC®(BBC) unit are in agreement with 
the previous work [16], in which the averaged a is prac¬ 
tically identical to the disordered B 13 C 2 , meanwhile the 
averaged a is slightly larger with respect to the exper¬ 
iments [ 8 , 42]. The ground state configuration, on the 
other hand, is more similar to the experiments for the 
parameter a but its angle a is slightly more deviated, as 
compared to those of the disordered configurations. As 
demonstrated by Isaev et al [44], the effects of lattice 
vibrations, i.e. zero-point motion and finite temperature 
effects, can play an important role in a system of light 
elements, e.g. elementary boron. In such a case, the in¬ 
clusion of vibrational effects corrects the underestimated 
equilibrium volume and gives excellent agreement with 
the experiment. We suggest that the underestimated lat¬ 
tice parameters for the disordered B 13 C 2 may be origi¬ 
nating from the neglect of lattice vibrations. Regarding 
the crystal symmetry, only the e(l, 4), the e(l- 6 ), and 
the e(l- 6 )+p(l- 6 ) configurations have, on average, the 
full rhombohedral symmetry with the R3m space group, 
while the e(l), and e(l-3) configurations are lacking in¬ 
version symmetry, thus yielding the R3m space group on 
average. We note that the existence of the R3m phase 
has been discussed by us in configurationally disordered 
B 4 C, stable at elevated temperatures [15], and more re¬ 
cently also by Yao et al. [45]. 

b. Bending and rotational effects in BBC chains 

After the ionic relaxation, we observe that the BBC 
chains in the disordered B 13 C 2 are bending, by which 
the chain-center boron atoms in the BBC chains are dis¬ 
placed to interstitial positions around their former posi¬ 
tions, resulting in empty space at the chain-center posi¬ 
tions. An illustration of the bended BBC chain is given 
by Eig. 5. We further investigate the bended BBC chain 
by 360° rotating the chain, and find that the chain-center 
boron atoms move toward and likely to form bonds with 
neighboring icosahedra bonded to the chain-end boron 
atoms. In addition, the chain-center boron atom is unfa¬ 
vorable to form bond with the icosahedral carbon atom 
bonded to the chain-end boron atom. These findings are 


in line with the work of Kwei et al. [46], where they 
demonstrated, using neutron diffraction technique, the 
existence of both the vacancies at the chain-center posi¬ 
tions and the interstitial atoms around the chain-center 
positions. They stated that this may originate from the 
presence of the rhombic B 4 or non-linear chains, such 
as BBB chain. However, we demonstrated in section III 
that the BBB chain is the high energy defect. As for the 
rhombic B 4 chain, it will be shown in the following section 
VLe that a model with the rhombic B 4 chain, proposed 
by Shirai et al. [19] is less favorable compared to our mod¬ 
els in the present work. We thus suggest that rather than 
the BBB and the B 4 chains, the vacancies at and the in¬ 
terstitial atoms around the chain-center positions could 
originate from bending of the BBC chains. Eurthermore, 
due to the chain bending, the distance between the two 
chain-end atoms in the BBC chain decreases and ranges 
between 2.45 and 2.65 A, which is in agreement with the 
value of 2.5 A reported by Kwei et al. [46]. Meanwhile 
the distance between the two chain-end boron atoms in 
the bending BBB chain and the rhombic B 4 chain in our 
case are 1.85 A and 2.85 A, respectively. 



FIG. 5: Illustration of a BBC chain before (grey spheres) and 
after (black spheres) ionic relaxations, and their neighboring 
icosahedra (white spheres). 

Apart from the BBC chains, we discuss about the CBC 
chains in the disordered B 13 C 2 . Unlike the bended BBC 
chains, the CBC chains are still linear after the relax¬ 
ation. The length of the linear CBC chains in the dis¬ 
ordered B 13 C 2 is slightly different from that of the CBC 
chain in the B 12 (CBC) unit ( 2.886 A) by less than 0.25%. 
Meanwhile the length of the CBC chain in the Bi 2 (CBC) 
unit increases by 0.52% approximately as compared to 
that of the CBC chain in the BiiC^(CBC) unit of B 4 C. 
This is in excellent agreement with the experimental re¬ 
sults existing in the literature [ 8 , 42, 47], in which the 
length of the CBC chains in B 13 C 2 increases approxi¬ 
mately from 0.21% to 0.34% with respect to that of the 
CBC chain in B 4 C. 








c. Configurational stability of the disordered B13C2 

We also determine the stability of substitutionally dis¬ 
ordered B 13 C 2 , with respect to the ordered ground state 
Bi 2 (CBC). The Gibbs free energy for disordered config¬ 
uration 7 at zero pressure is obtained from; 

-TS^, ( 2 ) 

where and are referred to the Gibbs free 

energy, the total energy, and the configurational en¬ 
tropy, respectively. T is the absolute temperature in 
Kelvin. Meanwhile the configurational entropy is cal¬ 
culated within a mean-field approximation, as given by; 

n 

S = -kBN'^Xi\n{xi), (3) 

i=l 

where is the Boltzmann constant. N and n are defined 
as the number of superatom sites in the supercell and 
the number of superatom types included in the super¬ 
cell, respectively. Xi refers to the concentration of type- 
i superatom. By investigating the rotational effects of 
the bended BBG chains in the disordered B 13 G 2 as men¬ 
tioned in section VI:b, we find that there are two ways of 
alignments that the bended chains can form bonds with 
the icosahedra and thus result in practically the same 
total energy. This bending-angle degeneracy provides a 
contribution to the configurational entropy. The extra 
configurational entropy from the bending and rotational 
degree of freedom of the BBG chains is thus included by 
assuming that the number of superatom types and the 
concentration for each superatom type, representing the 
BiiG®(BBG) units in Eq. 3 become twice and half, re¬ 
spectively. Fig. 6 illustrates the Gibbs free energy, mod¬ 
eled with the mean-field SA-SQS approach, for different 
substitutionally disordered B 13 G 2 plotted as a function 
of absolute temperature. As shown in TABLE III, the 
energy differences for the e(l), the e(l, 4), the e(l-3), 
and the e(l- 6 ) configurations, relative to the OS configu¬ 
ration, are approximately 0.25 eV/s.a., where the energy 
differences among them are small [on the order of 10 “^ 
eV/s.a.]. We find that, above 1546 K, the e(l- 6 ) configu¬ 
ration becomes stable, with respect to the ordered ground 
state. This is due to the higher configurational entropy 
for the e(l- 6 ) configuration, that lowers the Gibbs free 
energy at high temperature, more than those of the other 
types of disorder. The high-temperature phase of B 13 G 2 
can thus be represented by a mixture of the Bi 2 (GBG) 
and the BnG®(BBG) units, in which some chain-end 
carbon atoms G^^ (G^^) swap position with one of the 
neighboring equatorial boron atoms B®^, B®^, B®^ (B^^, 
B® 6 ) within the icosahedra bonded to the chain-end 
atoms. Even though the e(l- 6 )+p(l- 6 ) configuration has 
the highest configurational entropy, indicated by its high 
slope in Fig. 6 , among the others listed in TABLE III, 


its free energy is rather high even at high temperature. 
This corresponds to the results obtained from the dilute 
substitutional defects study in section III. 



FIG. 6: (Color online) Gibbs free energies for different sub¬ 
stitutionally disordered B 13 C 2 plotted as a function of abso¬ 
lute temperatures (at P = 0 GPa), relative to the Bi 2 (CBC) 
ground state (dashed line). 

So far only the e(l- 6 ) configuration, consisting of 
50% of the Bi 2 (GBG) and the BiiG®(BBG) units, i.e. 
(Bi 2 (CBG))o. 5 (BiiG®(BBG))o. 5 , is considered. We thus 
further investigate the e(l- 6 ) configuration, as generally 
denoted by (Bi 2 (GBC))i_^(BiiG^(BBC))^, by which 
we vary the concentration x of the BnG®(BBG) units. 
The concentration x, reflecting the degree of substitu¬ 
tional disorder in B13G2, is varied between 0 ( 100 % of 
Bi 2 (CBG)) and 1 ( 100 % of BiiG^(BBG)). Fig. 7 depicts 
a plot between the Gibbs free energy of the e(l- 6 ) config¬ 
uration and the concentration x at various absolute tem¬ 
peratures, where ^<x<\. It is worth noting again that 
the BnG®(BBG) units have 6 types of degeneracy, ac¬ 
cording to the substitution of the chain-end carbon atoms 
at one out of the six equatorial sites, and have 12 types 
of degeneracy if the rotational degree of freedom of the 
bended BBG chains is taken into account. In this study, 
the concentration of BnG®(BBG) is sampled discretely at 
x = 0, |, |, and 1. The results in Fig. 7 indicate 

that as the temperature increases, B13G2 is energetically 
favorable to disorder, where some of the Bi 2 (CBG) units 
will transform into the BnG®(BBC) units. According to 
our samplings, the e(l- 6 ) configuration with x = |, ^, 
and ^ will be energetically favorable when the temper¬ 
ature is higher than 1243 K, 1350 K, and 2350 K, re¬ 
spectively. Also, because of such a disorder, the electron 
deficiency in B13C2 will be gradually compensated until 
it eventually becomes a semiconductor, as x is approach¬ 
ing In our case, the hole density is reduced to 0.67, 
0.33, and 0 hole/s.a. for the e(l- 6 ) configuration with x 
= ^, ^, and respectively. As x is larger than the 
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stability of the disorderd B 13 C 2 tends to drastically de¬ 
crease, predicted by a relatively higher Gibbs free energy 
even at high temperatures. This is not only because of 
a constant increase of the less stable BnC® (BBC) units, 
but also a non-linear effect at (not shown). 



FIG. 7: (Color online) Plot between the Gibbs free energies 
of the e(l-6) configuration and the concentration x of the 
Bn(BBC) units at various absolute temperatures, relative 
to the Bi 2 (CBC) ground state (dashed line). 

Based on our results, a temperature as high as 2350 K 
is required in order to achieve a completely semiconduct¬ 
ing B 13 C 2 . We note that boron carbide can be manufac¬ 
tured, for instance, by hot-pressing, at typically above 
2000 K [1]. This high temperature corresponds to the 
study of Kuzenkova et al. [48] that the kinetics of recrys¬ 
tallization, reflecting atomic diffusion, in boron carbide 
starts at temperature above 2000 K due to its strong in¬ 
teratomic covalent bonds. Our results reveal, at 2000 K, 
the Gibbs free energy of the e(l- 6 ) configuration with x 
= ^ is higher than that with ^ | only by less than 1 

meV/atom, which is practically negligible. On the other 
hand, upon cooling, the atomic diffusion within boron- 
carbon compound might be limited/quenched, thus pre¬ 
venting it to reach the ordered metallic Bi 2 (CBG) state. 
We thus suggest that this could be the main cause for 
B 13 C 2 to behave as a semiconductor, i.e. due to the in¬ 
evitable substitutional disorder between boron and car¬ 
bon atoms. Our results corroborate the suggestions by 
Balakrishnarajan et al. [28], Schmechel [29], and Wer- 
heit [30], that semiconducting behavior in B 13 C 2 could 
arise from substitutional disorder. 


d. Estimation of the uncertainty in transition 
temperature predictions 

The transition temperatures, calculated in this work, 
involve several approximations that can give rise to the 


uncertainty in predicting transition temperatures. Con¬ 
sequently, we estimate in this section errors for the tran¬ 
sition temperatures, originating from different sources, 
and discuss how such errors affect the transition temper¬ 
atures. 

The first source of errors originates from the approx¬ 
imations in DFT itself. Since the exchange-correlation 
energy functionals are not known exactly, calculated re¬ 
sults, i.e. the transition temperature in our case, may 
either increase or decrease depending on used exchange- 
correlation functionals. Hautier et al. [49] recently mod¬ 
eled the computational errors within the DFT using GGA 
(+U for d-block metals) as the exchange-correlation func¬ 
tional. The errors were estimated by evaluating the 
reaction energies for ternary from binary oxides using 
DFT and compared them with experiments. The error, 
i.e. standard deviation, is thus approximated to be 24 
meV/atom. In our case, the error can be expected to be 
much smaller, since the bonding character and the local 
atomic environment are very similar for most atoms when 
comparing disordered boron carbides and the ordered 
ones. Meanwhile in the case of Hautier et a/., forming 
ternary oxides from different binary oxides does consid¬ 
erably change bonding characters as well as the atomic 
environment of most atoms. Instead, if we approximate 
that the DFT error in our case is around 24 meV per su¬ 
peratom, the structural unit observing distinct changes, 
we approximate the error bar for the order-disorder tran¬ 
sition temperature between the disordered and the or¬ 
dered B 13 C 2 , which are respectively the e(l- 6 ) and the 
ground-state Bi 2 (CBC), due to exchange-correlation un¬ 
certainties to be around 10%, or 150 K. 

Another source of errors is attributed to the use of the 
SA-SQS approach. Unlike normal atoms, each superatom 
has an internal structure, consisting of at least 15 atoms. 
Consequently, the commutative property is not preserved 
within the SA-SQS, i.e. BA % AB, in which A and B are 
superatoms of different types. To estimate the error, we 
consider four substitutionally disordered configurations 
of B4C, in which, for each configuration, the random sub¬ 
stitution of icosahedral carbon atoms C^ at all six polar 
sites coexists with the bipolar defects (B 10 C 2 +B 12 ) at 
high concentrations, as proposed to be the high temper¬ 
ature phase of B4C and being denoted by the notation 
BD+ 6 P in our previous work [15]. The four models of 
BD+ 6 P are obtained, using the same superatom basis 
and the same superatom types, as suggestes in Ref. [15], 
but they are allowed to interchange. Due to the com¬ 
mutation of superatom types, each BD+ 6 P has differ¬ 
ent internal configuration from the others, thus resulting 
in different total energies as well as different transition 
temperatures. We estimate the uncertainty by calculat¬ 
ing a mean value of transition temperature between the 
BD+ 6 P and the ground-state BiiC^(CBC) unit and its 
standard deviation. The mean transition temperature 
and the standard deviation are 1426 K and 152.45 K, 
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respectively. Therefore, the transition temperature of 
1426T305 K is obtained with 95% confidence intervals, 
i.e. the uncertainty in transition temperature caused by 
this approximation is < 20 %. 

The third error is due to the use of mean-field approx¬ 
imation for the configurational entropy. It is known that 
this approximation generally overestimates the order- 
disorder transition temperature by 20-30% due to ne¬ 
glecting short range ordering effects. An example is given 
by the work of Ailing [50], in which the order-disorder 
transition temperature of Tio. 5 Mgo. 5 N was calculated 
both with mean-field and Monte Carlo methods based on 
unified cluster expansion [51]. The mean-field transition 
temperature was found to be 1272 K, overestimating the 
more accurate Monte Carlo simulation that obtained the 
transition temperature of 950 K by 322 K or 30% approx¬ 
imately. Recently, Yao et al. [45] studied configurational 
phase transition in B 4 C by which they used Monte Carlo 
simulations. Based on their work, they proposed a pair of 
transitions. After the first transition at 717 K, they found 
that the icosahedral carbon atoms in BnC^ of the mono¬ 
clinic ground state B 4 C, represented by the BiiC^(CBC) 
units, are substitutionally disordering on the three polar- 
up (or down) sites, resulting in the RZm symmetry. This 
is in fact equivalent to the 3PU configuration in our pre¬ 
vious work of B 4 C [15]. Compared to ours, the transition 
temperature in our case is slightly higher by 150 K due 
to the use of mean-field approximation, i.e. overestima- 
tion with 20% approximately. For the higher transition, 
where icosahedral carbon atoms are disordering at all six 
polar sites (R3m), the difference in the transition tem¬ 
perature between the two works are larger, but further 
investigations are needed before this difference is referred 
only to the mean-field approximation. 

Perhaps the largest source of uncertainty in the present 
temperature prediction is the absence of explicit vibra¬ 
tional effects. We referred to the work of Isaev [44] in 
section VI:a that the effect of lattice vibrations can play 
an important role in a system of light elements. Garbul- 
sky et al. [52] demonstrated that the vibrational effect on 
order-disorder transition temperature can be significant, 
by which it tends to lower the transition temperature. 
Taking into account the vibrational effects is however be¬ 
yond the scope of this work. 

Taken together, our approach can give quantitative 
information about the properties of disordered phases, 
and a qualitative description of a series of order-disorder 
phase transitions, but not a quantitative description of 
their exact transition temperatures. It does however pro¬ 
vide the necessary starting point, in term of candidate 
structures, for future more elaborate investigations, in 
particular including vibrations. By taking into account 
the errors that could be arising from our approach, our 
transition temperatures are likely to be overestimated by 
several hundreds Kelvin. However, based on the discus¬ 
sion above, even though our transition temperatures are 


high, they are not going to exceed the melting tempera¬ 
ture of the material. As a result, the disordered B 13 C 2 
consisting of the Bi 2 (CBC) and the BiiC®(BBC) units 
can be a representation of a high temperature phase of 
B 13 C 2 and would be able to provide a solution to the 
electronic structure problem in B 13 C 2 . 



(a) 


(b) 


(c) 




FIG. 8: Three types of superatoms, used to construct struc¬ 
tural models for B13C2 proposed by Shirai et al. [ 19 ]: (a) 
BiiC^(CBC), (b) Bi2(CBC), and (c) 612(64) with a 4-atom 
rhombic chain. Light and dark spheres represent boron and 
carbon atoms, respectively. 


e. Comparison with models proposed in the 
literature 

Recently, Shirai et al. [19] proposed an alternative 
structural model for B 13 C 2 , consisting of 3 different types 
of structural units, i.e. BiiG^(CBG), Bi 2 (CBG), and 
312 ( 64 ) with a 4-atom rhombic ( 64 ) chain (see Fig. 8 ). 
We were thus inspired, with a use of SA-SQS, to exam¬ 
ine the model II in Ref. [19] in a 2x2x2 supercell with the 
same recipe (14.05 at.% C) for comparison. The model 
consists of 3 units of BiiC^(CBG), a unit of 612 ( 64 ), and 
4 units of Bi 2 (CBG). The total energy differences of the 
SA-SQS model II, relative to Bi 2 (CBG) and a chemical 
potential of diamond, are given in TABLE IV. Compared 
to the results in the literature, our results show a fairly 
good agreement. In fact, the total energies of the SA- 
SQS model II are somewhat lower due to the removal of 
cell constraints. However, having the icosahedral carbon 
atoms at the equatorial sites, bonded to the 4-atom rhom¬ 
bic chains (Model 11/(e)), in our case does not reduce the 
total energy as reported in the literature [19]. We under¬ 
line that the negative values of AE in TABLE IV are no 
proof of stability in thermodynamic sense. To examine 
their thermodynamic stability, one might need to con¬ 
sider also the energy difference with respect to B 4 C. The 
electronic DOS for the SA-SQS model H are illustrated 
by Eig. 9(a) and 9(b). Both of them have the same metal¬ 
lic character and the same hole density of 0.25 hole/s.a. 
as the model H in [19]. Even though the distances be¬ 
tween the valence band edge and the first gap state are 
different from that in the literature, the qualitatively cor¬ 
responding results for the model H are achieved with the 
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TABLE IV: Energy difference (AE) of Model II proposed in 
Ref. [19], relative to Bi 2 (CBC) plus a chemical potential of 
diamond, as a correction for one additional carbon atom. As 
given in Ref. [19], type (p) indicates configurations, in which 
all icosahedral carbon atoms are located at the polar site, 
meanwhile (e) indicates configurations, in which at least one 
of the icosahedral carbon atoms is located at the equatorial 
site. 


Gonfiguration/(type) 

AE (eV/atom) 

Model II/(p) 

-0.004 “ 


0.007 (Ild) 


0.011 (lie) ’’ 

Model II/(e) 

-0.004 “ 


0.004 (Ilf) 


0.001 (Ilg) *’ 


“This work (GGA-PAW without cell constraint) 

^Ref. [19] - Shirai at al. (LDA-PP with cell constraints) 

SA-SQS approach, strengthening the reliability of the ap¬ 
proach. 



EIG. 9: Density of states of (a) Model II/(p), (b) Model 
II/(e), and (c) Rh-Bi 3 C 2 , obtained by using the GGA-PBE96 
for the exchange-correlation functional. The Eermi level is lo¬ 
cated at 0 eV. 

We then construct another SA-SQS model of B 13 C 2 , 
based the three types of superatom in Fig. 8 , in a 4x4x3 
supercell. The cell consists of 48 superatoms (735 atoms): 
32 units of BiiC^(CBC), 15 units of Bi 2 (B 4 ), and a 
unit of Bi 2 (CBC). This is also inspired by the sugges¬ 
tion of Shirai et al. [19], but in contrast to that work, our 
larger supercell contains exact stoichometry of B 13 C 2 . As 
shown in Fig. 9(c), this SA-SQS model of B 13 C 2 , namely 
Rh-Bi 3 C 2 , is a semiconductor with a GGA-PBE96-based 
band-gap of 1.42 eV. The formation energy, calculated 
from a-boron and diamond is -0.049 eV/atom. This is 
higher than that of (Bi 2 (CBC))o. 5 (BiiC®(BBC))o .5 of 
the e(l- 6 ) configuration by 0.027 eV/atom at 0 K. We 


TABLE V: Formation energy (AE-^^’^"^) of different 15-atom 
B 14 G units, with respect to o-boron and diamond. 


Structural unit 

(eV/atom) 

Bi 2 (BBC) 

0.010 

BiiCP(BBB) 

0.067 

BiiC'^(BBB) 

0.078 

Bi 2 (BCB) 

0.211 


note that to determine the configurational entropy, in 
this case, we take into account the possibility of hav¬ 
ing ( 1 ) the substitution of icosahedral carbon atoms at 
different polar sites and ( 2 ) three fold rotational degen¬ 
eracy of the 4-atom rhombic chains to obtain an upper 
bound of the configurational entropy. Nevertheless, Rh- 
B 13 C 2 is not favorable with respect to the e(l- 6 ) config¬ 
uration at any realistic temperature below the extreme 
17000 K. The total energy of Rh-Bi 3 C 2 can be lowered, 
as shown in Ref. [19], by increasing a concentration of the 
Bi 2 (CBC) ground state units. However, one can expect 
that by doing so the compound would be more metal¬ 
lic. Consequently, based on our results, we propose that 
at 13.33 at.% C, the atomic configuration of boron car¬ 
bide is dominated by a mixture of the Bi 2 (CBC) and the 
BiiC®(BBC) units, although other structural units, like 
BiiC^(BBC), BiiC^(CBC) and Bi 2 (B 4 ) can be present 
at lower concentrations. 


VII. IMPLICATIONS ON UNDERSTANDING 
OTHER BORON CARBIDES 

Based on our model proposed in this work, as well 
as our previous study of B 4 C [15], we can suggest the 
following scenario for structural evolution as the stoi¬ 
chiometry is deviating from the two "ideal” B 13 C 2 and 
B 4 C compositions: going from B 13 C 2 toward the carbon- 
rich limit (~20 at.% C, or B 4 C), the Bi 2 (CBC) and the 
BiiC®(BBC) units would gradually be replaced by the 
stable, semiconducting BiiC^(CBC) units with config¬ 
urational disorder of carbon atoms on the icosahedral 
polar sites [15], until the latter dominates the boron car¬ 
bide system at the carbon-rich limit as suggested in the 
literature [12-15]. On the other hand, if the carbon con¬ 
centration becomes lower than 13.33 at.% C toward the 
boron-rich limit (~8 at.% C), an equal fraction of car¬ 
bon atoms in the BuC® icosahedra and the (CBC) chains 
could be substituted by boron atoms, thus yielding the 
Bi 2 (BBC) units, which we find to have the lowest forma¬ 
tion energy ( 0.01 eV/atom), among the other B 14 C units 
(see TABLE V), with respect to a-boron and diamond. 
As shown in section IV, the Bi 2 (BBC) unit is a semi¬ 
conductor. Thus, based on this structural model, boron 
carbide would maintain semiconducting character within 
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its whole single-phase region due to specific types of low 
energy configurational disorder stabilized by entropy as 
shown in detail for B 13 C 2 in the present work. However, 
apart from the stoichiometries B4C and B13C2, details 
of the stability and the properties of other boron carbide 
compositions deserve further investigation. 

VIII. CONCLUSIONS 

We investigate the properties of B 13 C 2 using first- 
principles calculations, in which substitutionally disor¬ 
dered configurations are modelled within the SA-SQS 
scheme. The calculations of dilute stoichiometric defects 
reveal, in the matrix of the presumed ground state B 13 C 2 
represented by the Bi 2 (CBC) units, the BiiC®(BBC) 
unit have the lowest formation energy. We thus predict 
that as the temperature increases, the low temperature- 
ordered B 13 C 2 , comprising only of the Bi 2 (CBC) units, 
becomes thermodynamically unfavorable with respect to 
substitutional disorder, in which some of the Bi 2 (CBC) 
units turn into the BnC®(BBC) units, leading to sub¬ 
stitutionally disordered B 13 C 2 . The splitting of valence 
states in the BnC®(BBC) units then compensates par¬ 
tially or fully the electron deficiency in the Bi 2 (CBC) 
units, depending on their relative concentrations. Con¬ 
sequently, as the concentration of the BnC® (BBC) units 
increases with temperature under equilibrium condition, 
approaching 50%, the electron deficiency in B 13 C 2 is 
gradually compensated until the completely semicon¬ 
ducting character is achieved. Also the calculated band- 
gap for disordered B 13 C 2 is in fairly good agreement with 
the experiment. It is possible that the metallic ordered 
B 13 C 2 has not been experimentally observed because of 
the limited atomic diffusion in boron carbide even at ele¬ 
vated temperature, thus freezing in the high temperature 
disordered configurations of B 13 C 2 also at intermediate 
temperature. 

This mixture of the Bi 2 (CBC) and the BiiC®(BBC) 
units is found to have lower total energy and Gibbs 
free energy than a model based on previously suggested 
BiiC^(CBC), Bi 2 (B 4 ), and Bi 2 (CBC) units [19]. The 
model for structural disorder proposed in the present 
work can explain experimental finding of semiconduct¬ 
ing character of B 13 C 2 . 
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